library(raster)
library(rgeos)
library(rgdal)
library(sp)
library(sf)
library(velox)

load("data/grid_data.Rdata")

######0 AD

# Read raster
pop0rast <- raster("data/Hydepop/hyde31_final/0ad_pop/popc_0AD.asc")
# Set NA pixels to 0
pop0rast[is.na(pop0rast)] <- 0
# Convert to velox object
pop0velox <- velox(pop0rast)
# Extract raster to polygon
cntrypopsum <- pop0velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e_priogrid@data$gid,cntrypopsum)
rm(pop0rast, pop0velox)


######100 AD

# Read raster
pop100rast <- raster("data/Hydepop/hyde31_final/100ad_pop/popc_100AD.asc")
# Set NA pixels to 0
pop100rast[is.na(pop100rast)] <- 0
# Convert to velox object
pop100velox <- velox(pop100rast)
# Extract raster to polygon
cntrypopsum <- pop100velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop100rast, pop100velox)

######200 AD

# Read raster
pop200rast <- raster("data/Hydepop/hyde31_final/200ad_pop/popc_200AD.asc")
# Set NA pixels to 0
pop200rast[is.na(pop200rast)] <- 0
# Convert to velox object
pop200velox <- velox(pop200rast)
# Extract raster to polygon
cntrypopsum <- pop200velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop200rast, pop200velox)

######300 AD

# Read raster
pop300rast <- raster("data/Hydepop/hyde31_final/300ad_pop/popc_300AD.asc")
# Set NA pixels to 0
pop300rast[is.na(pop300rast)] <- 0
# Convert to velox object
pop300velox <- velox(pop300rast)
# Extract raster to polygon
cntrypopsum <- pop300velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop300rast, pop300velox)

######400 AD

# Read raster
pop400rast <- raster("data/Hydepop/hyde31_final/400ad_pop/popc_400AD.asc")
# Set NA pixels to 0
pop400rast[is.na(pop400rast)] <- 0
# Convert to velox object
pop400velox <- velox(pop400rast)
# Extract raster to polygon
cntrypopsum <- pop400velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop400rast, pop400velox)

######500 AD

# Read raster
pop500rast <- raster("data/Hydepop/hyde31_final/500ad_pop/popc_500AD.asc")
# Set NA pixels to 0
pop500rast[is.na(pop500rast)] <- 0
# Convert to velox object
pop500velox <- velox(pop500rast)
# Extract raster to polygon
cntrypopsum <- pop500velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop500rast, pop500velox)


######600 AD

# Read raster
pop600rast <- raster("data/Hydepop/hyde31_final/600ad_pop/popc_600AD.asc")
# Set NA pixels to 0
pop600rast[is.na(pop600rast)] <- 0
# Convert to velox object
pop600velox <- velox(pop600rast)
# Extract raster to polygon
cntrypopsum <- pop600velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop600rast, pop600velox)

######700 AD

# Read raster
pop700rast <- raster("data/Hydepop/hyde31_final/700ad_pop/popc_700AD.asc")
# Set NA pixels to 0
pop700rast[is.na(pop700rast)] <- 0
# Convert to velox object
pop700velox <- velox(pop700rast)
# Extract raster to polygon
cntrypopsum <- pop700velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop700rast, pop700velox)

######800 AD

# Read raster
pop800rast <- raster("data/Hydepop/hyde31_final/800ad_pop/popc_800AD.asc")
# Set NA pixels to 0
pop800rast[is.na(pop800rast)] <- 0
# Convert to velox object
pop800velox <- velox(pop800rast)
# Extract raster to polygon
cntrypopsum <- pop800velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop800rast, pop800velox)


######900 AD

# Read raster
pop900rast <- raster("data/Hydepop/hyde31_final/900ad_pop/popc_900AD.asc")
# Set NA pixels to 0
pop900rast[is.na(pop900rast)] <- 0
# Convert to velox object
pop900velox <- velox(pop900rast)
# Extract raster to polygon
cntrypopsum <- pop900velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop900rast, pop900velox)

######1000 AD

# Read raster
pop1000rast <- raster("data/Hydepop/hyde31_final/1000ad_pop/popc_1000AD.asc")
# Set NA pixels to 0
pop1000rast[is.na(pop1000rast)] <- 0
# Convert to velox object
pop1000velox <- velox(pop1000rast)
# Extract raster to polygon
cntrypopsum <- pop1000velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop1000rast, pop1000velox)

######1100 AD

# Read raster
pop1100rast <- raster("data/Hydepop/hyde31_final/1100ad_pop/popc_1100AD.asc")
# Set NA pixels to 0
pop1100rast[is.na(pop1100rast)] <- 0
# Convert to velox object
pop1100velox <- velox(pop1100rast)
# Extract raster to polygon
cntrypopsum <- pop1100velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop1100rast, pop1100velox)


######1200 AD

# Read raster
pop1200rast <- raster("data/Hydepop/hyde31_final/1200ad_pop/popc_1200AD.asc")
# Set NA pixels to 0
pop1200rast[is.na(pop1200rast)] <- 0
# Convert to velox object
pop1200velox <- velox(pop1200rast)
# Extract raster to polygon
cntrypopsum <- pop1200velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop1200rast, pop1200velox)



######1300 AD

# Read raster
pop1300rast <- raster("data/Hydepop/hyde31_final/1300ad_pop/popc_1300AD.asc")
# Set NA pixels to 0
pop1300rast[is.na(pop1300rast)] <- 0
# Convert to velox object
pop1300velox <- velox(pop1300rast)
# Extract raster to polygon
cntrypopsum <- pop1300velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop1300rast, pop1300velox)




######1400 AD

# Read raster
pop1400rast <- raster("data/Hydepop/hyde31_final/1400ad_pop/popc_1400AD.asc")
# Set NA pixels to 0
pop1400rast[is.na(pop1400rast)] <- 0
# Convert to velox object
pop1400velox <- velox(pop1400rast)
# Extract raster to polygon
cntrypopsum <- pop1400velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop1400rast, pop1400velox)






######1500 AD

# Read raster
pop1500rast <- raster("data/Hydepop/hyde31_final/1500ad_pop/popc_1500AD.asc")
# Set NA pixels to 0
pop1500rast[is.na(pop1500rast)] <- 0
# Convert to velox object
pop1500velox <- velox(pop1500rast)
# Extract raster to polygon
cntrypopsum <- pop1500velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop1500rast, pop1500velox)





######1600 AD

# Read raster
pop1600rast <- raster("data/Hydepop/hyde31_final/1600ad_pop/popc_1600AD.asc")
# Set NA pixels to 0
pop1600rast[is.na(pop1600rast)] <- 0
# Convert to velox object
pop1600velox <- velox(pop1600rast)
# Extract raster to polygon
cntrypopsum <- pop1600velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop1600rast, pop1600velox)





######1700 AD

# Read raster
pop1700rast <- raster("data/Hydepop/hyde31_final/1700ad_pop/popc_1700AD.asc")
# Set NA pixels to 0
pop1700rast[is.na(pop1700rast)] <- 0
# Convert to velox object
pop1700velox <- velox(pop1700rast)
# Extract raster to polygon
cntrypopsum <- pop1700velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop1700rast, pop1700velox)




######1800 AD

# Read raster
pop1800rast <- raster("data/Hydepop/hyde31_final/1800ad_pop/popc_1800AD.asc")
# Set NA pixels to 0
pop1800rast[is.na(pop1800rast)] <- 0
# Convert to velox object
pop1800velox <- velox(pop1800rast)
# Extract raster to polygon
cntrypopsum <- pop1800velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop1800rast, pop1800velox)

######1900 AD

# Read raster
pop1900rast <- raster("data/Hydepop/hyde31_final/1900ad_pop/popc_1900AD.asc")
# Set NA pixels to 0
pop1900rast[is.na(pop1900rast)] <- 0
# Convert to velox object
pop1900velox <- velox(pop1900rast)
# Extract raster to polygon
cntrypopsum <- pop1900velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop1900rast, pop1900velox)

######2000 AD
# Read raster
pop2000rast <- raster("data/Hydepop/hyde31_final/2000ad_pop/popc_2000AD.asc")
# Set NA pixels to 0
pop2000rast[is.na(pop2000rast)] <- 0
# Convert to velox object
pop2000velox <- velox(pop2000rast)
# Extract raster to polygon
cntrypopsum <- pop2000velox$extract(sp = e_priogrid, fun=sum, df = TRUE, small = TRUE)
e1 <- cbind(e1,cntrypopsum$out)
rm(pop2000rast, pop2000velox)

# names(e1)

####merging on natharbdist
uniquesetVdem$natharbdist_km <- uniquesetVdem$two_port_all_pred1/1000

uniquesetVdem$natharbdist_km <- uniquesetVdem$natharbdist_km/100

d <- subset(uniquesetVdem, year==2005, select=c("gid", "natharbdist_km"))

colnames(e1)[1] <- "gid"

e1 <- merge(e1, d, by=c("gid"), all.x=T)
rm(uniquesetVdem)


sum(e1[,23])
####REGRESSIONS

m1 <- lm(e1[,3] ~ e1$natharbdist_km)
m2 <- lm(e1[,4] ~ e1$natharbdist_km)
m3 <- lm(e1[,5] ~ e1$natharbdist_km)
m4 <- lm(e1[,6] ~ e1$natharbdist_km)
m5 <- lm(e1[,7] ~ e1$natharbdist_km)
m6 <- lm(e1[,8] ~ e1$natharbdist_km)
m7 <- lm(e1[,9] ~ e1$natharbdist_km)
m8 <- lm(e1[,10] ~ e1$natharbdist_km)
m9 <- lm(e1[,11] ~ e1$natharbdist_km)
m10 <- lm(e1[,12] ~ e1$natharbdist_km)
m11 <- lm(e1[,13] ~ e1$natharbdist_km)
m12 <- lm(e1[,14] ~ e1$natharbdist_km)
m13 <- lm(e1[,15] ~ e1$natharbdist_km)
m14 <- lm(e1[,16] ~ e1$natharbdist_km)
m15 <- lm(e1[,17] ~ e1$natharbdist_km)
m16 <- lm(e1[,18] ~ e1$natharbdist_km)
m17 <- lm(e1[,19] ~ e1$natharbdist_km)
m18 <- lm(e1[,20] ~ e1$natharbdist_km)
m19 <- lm(e1[,21] ~ e1$natharbdist_km)
m20 <- lm(e1[,22] ~ e1$natharbdist_km)
m21 <- lm(e1[,23] ~ e1$natharbdist_km)


coefs <- rbind(coef(m1)[2],
               coef(m2)[2],
               coef(m3)[2],
               coef(m4)[2],
               coef(m5)[2],
               coef(m6)[2],
               coef(m7)[2],
               coef(m8)[2],
               coef(m9)[2],
               coef(m10)[2],
               coef(m11)[2],
               coef(m12)[2],
               coef(m13)[2],
               coef(m14)[2],
               coef(m15)[2],
               coef(m16)[2],
               coef(m17)[2],
               coef(m18)[2],
               coef(m19)[2],
               coef(m20)[2],
               coef(m21)[2])


ses<- rbind(coef(summary(m1))[2, "Std. Error"],
            coef(summary(m2))[2, "Std. Error"],
            coef(summary(m3))[2, "Std. Error"],
            coef(summary(m4))[2, "Std. Error"],
            coef(summary(m5))[2, "Std. Error"],
            coef(summary(m6))[2, "Std. Error"],
            coef(summary(m7))[2, "Std. Error"],
            coef(summary(m8))[2, "Std. Error"],
            coef(summary(m9))[2, "Std. Error"],
            coef(summary(m10))[2, "Std. Error"],
            coef(summary(m11))[2, "Std. Error"],
            coef(summary(m12))[2, "Std. Error"], 
            coef(summary(m13))[2, "Std. Error"],
            coef(summary(m14))[2, "Std. Error"],
            coef(summary(m15))[2, "Std. Error"],
            coef(summary(m16))[2, "Std. Error"],
            coef(summary(m17))[2, "Std. Error"],
            coef(summary(m18))[2, "Std. Error"],
            coef(summary(m19))[2, "Std. Error"],
            coef(summary(m20))[2, "Std. Error"],
            coef(summary(m21))[2, "Std. Error"])

years <- as.matrix(seq(from=0, to=2000, by=100))


tiff("output/figure_7_2.tiff")
plot(NA, xlim = c(1, 21) , ylim = c(-7000, 2000), xlab = "", ylab = "", yaxt = "n", xaxt="n")
axis(1, 1:21,years, las = 2)
axis(2, seq(from=2000, to=-7000, by=-200), las = 2)
axis(4, -2500, labels="Coefficient", tck=0)
axis(3, 10.5, labels="Year", tck=0)
# We'll add a vertical line for zero:
abline(v = 0, col = "gray")
points( 1:21,coefs[,1], pch = 23, col = "grey", bg = "grey")
for(i in 1:21){
  segments(i, (coefs[i,1] - qnorm(0.975) * ses[i,1]), i, (coefs[i,1] + qnorm(0.975) *ses[i,1]), col = "black", lwd = 2)
}
dev.off()

